Anomalous diffusion in correlated continuous time random walks 



O 

(N 

q 

>—> 

in 

(N 



I 

C 

o 
o 



ON 
O 



X 



Vincent Tejedor 1 and Ralf Metzler 1 

Physics Department T30g, Technical University of Munich, 85747 G arching, Germany 

We demonstrate that continuous time random walks in which successive waiting times are 
correlated by Gaussian statistics lead to anomalous diffusion with mean squared displacement 
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Long-ranged correlations of the waiting times with power-law exponent a (0 < a < 2) 



give rise to subdiffusion of the form (r 2 (i)} ~ f ex /( 1 + a )_ j n contrast correlations in the jump lengths 
are shown to produce superdiffusion. We show that in both cases weak ergodicity breaking occurs. 
Our results are in excellent agreement with simulations. 



PACS numbers: 02.50.Ey,05.40.- 



I. INTRODUCTION 



The continuous time random walk (CTRW) theory 
was introduced more than forty years ago [l| to ex- 
tend regular random walks on lattices to a continuous 
time variable. We characterise each jump in a CTRW 
process by a waiting time r elapsing after the previ- 
ous jump, and a variable jump length £. Each r and 
£ are independent random variables, identically distrib- 
uted according to the probability densities ip(r) and A(£). 
For a power-law form ip(r) ~ T - x -& with < j3 < 1 
the characteristic waiting time J Q Tip(r)dT diverges and 
the resulting CTRW is subdiffusive, (r 2 {t)} ~ t p . For 
A(£) ~ |£| _1 1 with < 7 < 2 the jump length variance 
£ 2 A(£)d£ diverges and we obtain a superdiffusive 
Levy flight. Spatiotemporal coupling of jump length and 
waiting time leads to Levy walks with finite (r 2 (t)) ~ t 13 
with 1< (3 < 2 |, The CTRW model was cham- 
pioned in the seminal work on charge carrier transport in 
amorphous semiconductors [3]. CTRW theory has also 
been successfully applied in subsurface tracer dispersion 
[H, tick-tick dynamics in financial markets [f|, cardiolo- 
gical rhythms @, electron transfer Q, noise in plasma 
devices [811, dispersion in turbulent systems [9fl, search 
models [10|], or in models of gene regulation among 
many others [l2j]. CTRW models are closely related to 
the fractional Fokker-Planck equation [12, [Hj] ■ 

The independence of the waiting times and jump 
lengths giving rise to a renewal process is not always jus- 
tified. As soon as the random walker has some form of 
memory, even a short one, the variables become non- 
independent. Examples are found in financial market 
dynamics single trajectories in which there is a 

directional memory [l5j |. or in astrophysics [Hj]. An 
important application are search processes and human 
motion patterns in which memory and conscience will 
likely lead to a non-renewal situation. A general math- 
ematical framework was developed for non-independent 
CTRWs [13], however, it is quite cumbersome to apply 
and the cases solved so far only lead back to normal dif- 
fusion. An approach to coupling of waiting times based 
on the Langevin equation formulation of CTRW pro- 
cesses was recently introduced fl8| . Some special cases 
were explored that bridge between CTRW and fractional 



Brownian motion [l9| . Here we introduce a simple way to 
establish correlations in CTRW processes. Correlations 
between successive waiting times are shown to give rise to 
subdiffusion even when they are Gaussian, while correla- 
tions between jump lengths produce superdiffusion. We 
also consider long-ranged correlations and discuss anti- 
correlation effects. Our scaling arguments are in excel- 
lent accord with simulations. 

After introducing the general framework we demon- 
strate how correlated waiting times lead to subdiffusion 
and weak ergodicity breaking. We then consider correla- 
tions in the jump length, and proceed to analyse anticor- 
relation effects. Finally we provide some details on the 
derivations and draw our conclusions. 



II. GENERAL FRAMEWORK 

Let us briefly review the general framework for correl- 
ated CTRWs from Ref. 17] for correlated waiting times. 
Then the waiting time ip n for a given step n is conditioned 
by the previous waiting time ip n -i, as quantified by the 
joint probability P(ip n ,ip n _i) = P(ij! n \ip n _i)P(iji n _i). 
Here P(ip n -i) is the probability to have the waiting time 
ij) n -l in the previous step, and P(ip n \ip n -i) is the con- 
ditional probability to find a waiting time i\) n for given 
ipn-i- The normalisation conditions are 



OO />OG />OG 

P(lp n )dl(j n = 1, / / P(V>„,V>n-l)#n#n-l = 1) 
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P{lp n \i>n-{)dl]} n = 1. 



By recurrence, we obtain the joint probability 



(1) 



P (ipn, Ipn 



^o)=n P (^l^-lWo), (2) 



fc=l 



demonstrating that the waiting time ip n in fact depends 
on all previous waiting times. Note that in the decoupled 
case P(ij) k \ip k _i) = f{il>k)g(i>k-i) we get back to a regu- 
lar renewal CTRW. 

The marginal probability of ip n is defined as P(ip n ) = 
J °° P(ipn\ipn-i)P(iJn-i)diJn-i- According to Eq. ©, 
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this leads to an n-fold integration over the product on 
its right hand side. In the general case, it is quite hard 
to compute this quantity, and this is why in previous lit- 
erature only normal diffusion was treated in this frame- 
work. 



where we assumed r oc a. At long times the Gaussian 
waiting time correlation results in the subdiffusion law 

El 



(v 2 (t)) = (Sr 2 )(n(t))~Kt 
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(8) 



III. RANDOM WALK OF WAITING TIMES 

Instead of constructing the process from definitions ([T]) 
and @ we start from a different angle. Namely we build 
the value of waiting time ip n from the previous waiting 
time plus a random deviation Sip n , 



Ipn = i>n-l + Slp n . 



(3) 



The sequence of waiting times can therefore be viewed as 
a random walk in the space of waiting times. Similarly 
we will proceed with coupled jump lengths below. The 
waiting time ip n can therefore be expressed through ip n = 
Si=o where we assumed that ipo = Stpo, without loss 
of generality. A reflecting boundary condition at ip = 
ensures that the waiting times are positive. If the random 
variations Sip n are normally distributed, we obtain the 
following conditional probability: 



P(lp n \lpn- 



l) = 



exp 



exp 



Or, 



4cr 2 



4ct 2 



(4) 



The mean squared displacement (MSD) for this process 
grows with the number of steps as ((Aip(n)) 2 ) ~ 2er 2 n for 
large n. To proceed, we compute the probability P(t, ri) 
to have made n steps up to time t. In Laplace space [20[ 



P(s,n) 



P(t,n)e- St dt = L^Mfl TT 
" s to 



(5) 



Here and in the following we use the convention that 
the Laplace transform of a function f(t) is expressed 
by explicit dependence on the Laplace variable, that is, 
-^{/(*)} — fi s )- After some calculations we arrive at 
(see Appendix |A"| 



P(s,n) 



^u/ n(n+1K2n+1) , 



-ty I ^/(" + l)(" + 2)(2n + 3) a 



(6) 



We obtain the Laplace transform of the mean number 
of steps by summation, (n(s)) = J2^L nP{s,n). With 
the approximations detailed in Appendix [XJ we find in 
leading order around s = (corresponding to long times) , 
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where (Sr 2 ) is the jump length variance of the process 
and K = (Sr 2 )/a 2 ^ 3 the generalised diffusion constant. 
Eq. (JSJ) is one of the main result of this work. We stress 
again that this subdiffusion emerges from a random pro- 
cess that has a finite waiting time in each step. How- 
ever, as time proceeds this waiting time slowly diverges 
[(Aip(n)) ~ n 1 / 2 ], due to the diffusive coupling of waiting 
times. 

If we take a sharp correlation of waiting times with a = 
0, Eq. © leads to P(tp n ) = 5{ip n - tp ) + 5{ip n + tp ). At 
each step the waiting time is the same. The mean wait- 
ing time is (ip n ) = / n °° ip[5(ip- ip ) + 8(ip + ipo)] dip = ip 
such that we find classical Brownian motion with diffu- 
sion coefficient (5r 2 )/ipo, as it should be. 

Consider now the case of a stable distribution with 
index a for Sip. For 1 < a < 2 the first moment {Sip) = r 
of the random walk in ip space is still finite. We then 
follow similar steps as outlined above, obtaining 

(n(t)) - (t/a) a/{a+1) => (r 2 (t)) ~ Kt a ^ a+1 \ (9) 

The case < a < 1 is somewhat more involved. 
We argue that for any stable distribution, we have 

E"=i Mt^n-^+^^Sipin-^+^^t), where ~ is a scal- 
ing equality of distributions. This leads to the relation 
t(n)/cr~n ( ' a + 1 ^ a , such that Eq. © still holds for any 
< a < 2. We note that result © was retrieved from a 
Langevin equation approach in Rcf. 18]. 

For a stable distribution the jumps between successive 
waiting times become increasingly larger when a is de- 
creased, effecting even slower diffusion (r 2 (i)) ~ 
In the limit a — 2 we are back to Gaussian diffusion 
in ip space and a 2/3 subdiffusion in position space, 
(r 2 (i)) ~ t 2 l z . In Fig. [T]we demonstrate excellent agree- 
ment between our analytical findings and simulations res- 
ults of random processes performed with stable correla- 
tions between successive waiting times. 

The characteristic function of the stable variable Sip is 
'Psipil) — exp (—\cq\ a ), see Refs. [2l|. While at long times 
the predicted subdiffusive behaviour is attained, we ob- 
serve a transient regime of normal diffusion, (r 2 (t)) ~ t, 
when the scale factor c increases (see Fig. [2]) . At short 
times, we may neglect the probability that the random 
walker makes more than one step. The probability to 
make the first step corresponds to the cumulative func- 
tion of the waiting time distribution Fs^(t), 

2 f°° sin(to) , . . 2r(l + 1/a) , s 
F H {t) = - \ —±*lfai,(q)dq i '—Lt. (10) 

7T J Q q 7TC 

Thus, the initial linear slope in the ensemble average is 
due to the linearity of the cumulative function for short 
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Figure 1: (r 2 (i)} for a waiting time correlated 3D Gaussian 
walk. The Sip follow an a-stable law with scale factor c = 
1; q decreases from top to bottom. Simulations ( — ) and 
power-laws (• ■ • ) with fitted exponents 0.35, 0.50, 0.60, 0.66. 
Theoretical values a/ (a + 1): 0.33, 0.50, 0.60, 0.66. 
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Figure 2: Same as Fig. \T\ with a — 0.5 and various scale 
factors c (c increases from top to bottom). Simulations ( — 
) and power-laws (•) with slope 0.33 (c = 1, 10, 100) and 1 
(c = 1.000, 10.000). Note the turnover to slope a/(a + 1) at 
t ~ c. Compare text. 



waiting times. This effect vanishes as soon as the probab- 
ility increases that the walker makes two or more steps, 
and converges to the predicted long time behaviour (|9|). 



A. Weak ergodicity breaking 

When dealing with the time series of a single particle 
trajectory of length T, instead of the ensemble averaged 
MSD (r 2 (t)) one calculates the time averaged MSD 



-'(A,T) = 



1 



T — A 



T-A 



r(i + A) - r(t) dt, (11) 



relating two positions separated by the lag time A. 
Ensemble averaging Eq. (fTTj) the square brackets be- 




Figure 3: Time-averaged MSD as function of A for a CTRW 
with correlated waiting time (T = 2000). The plots show 
a scatter between individual trajectories that becomes more 
pronounced for decreasing a. The data approximately scale as 
the ensemble average, Eq. (| 13f) . For very small a we observe a 
plateau due to the occurrence of extremely long waiting times 
of the order of T. 



is the (finite) variance of jump lengths, and (nt,t+A) 
counts the average number of jumps in the time inter- 
val [t,t + A]. For normal diffusion (n(t)) = t/r and 
therefore <5 2 (A,T) = 2dKA behaves exactly as the en- 
semble average such that the lag time A is exchangeable 
with the process time t of the ensemble average. In con- 
trast for a subdiffusive renewal CTRW with waiting time 
density ip(t) ~ t^ 1 ^" (0 < a < 1) it turns out that 
5 2 (A,T) ~ KA/T 1 ~ a , i.e., we observe a so-called weak 
ergodicity breaking for A -C T [H, [23| ■ Do we observe 
similar behaviour for our coupled CTRW? With Eq. ([7]) 
and the relation (nt.t+A) = (™o,t+A> - (no.t), 



?(t + A) -r(t)Y 



2K 



(t + A) a /( 1+Q ) — t a /( 1 + a ) 

(12) 

with K = (<5r 2 )/(2dcr Q /( 1+Q )), and therefore 



<5 2 (A,T)) - 2dKA/T 



a/(l+a) 



(13) 



come ([r(i + A) - r(t)Y 



((5r 2 )(n tjt+A ) where (5r 2 



Thus our non-renewal, coupled CTRW also exhibits a 
weak ergodicity breaking, even in the limit a — 2 of 
Gaussian waiting time coupling. Simulations of this pro- 
cess indeed confirm the predicted scaling of the time av- 
eraged MSD with both lag time A and measurement time 
T. In Fig. [3] we demonstrate that individual trajectories 
show significant scatter in their amplitude, as observed 
for renewal CTRW subdiffusion [H. 



IV. CORRELATED JUMP LENGTHS 

Let us now consider a random walk with constant wait- 
ing time ■0 = 1 but correlated jump lengths. If r(t) is the 
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Figure 4: 3D random walk with correlated jump lengths. We 
chose x(0) = 0. 8x(t) is normally distributed with a = \f2. 
Simulations (o) and theory ( — ), Eq. (115[l , 



position of the walker we define x(i) = r(t) — r(t — I) 
(x(0) = 0). The jump length is now assumed to diffuse 
in x space, with increments Sx(t) = x(i) — x(< — I ) that 
are normally distributed with mean and variance a. We 
find 



P 



1 



{ira 2 t) 3 / 2 



exp 



x 2 + y 2 + z 2 



a 2 t 



(14) 



-a 2 t. For this process we have 



with variance (x 2 (t)) - ., 
r(t) = r(0) + Yli=i x (*) an( l obtain (see Appendix |A1 



t{t + \){2t + 1)(T 2 



(15) 



Thus, when we assume a diffusion in the space of jump 
lengths we obtain a Richardson-type t 3 scaling. Fig. @] 
demonstrates excellent agreement of Eq. (jf 51) with the 
simulations result. For a fixed jump length (a — 0) all 
steps have the same length and direction, and the corres- 
ponding random walk is ballistic, (r 2 (i)) = x 2 (0)£ 2 . 

Classical diffusion cannot be retrieved with this mech- 
anism, as this would require directional randomness for 
jumping left/right. In fact we obtain for the jump cor- 
relations that (x(i) • x(t 
correlations become 



t)) = \o 2 t. The position 



(r(t).r(t+r)) 



t(t + l)(2f+ \)a 2 3a 2 t(t + iy, 



(16) 



If the increments Sx are distributed according to a 
stable distribution of exponent a and scale factor c, 
each coordinate of r(t) is distributed according to a 
stable distribution of exponent a and scale factor c(t) — 

C : 



^Si=i *°J c — ct^ 1+a ^ Q . The mean squared displace- 
ment diverges for < a < 2 but we observe the super- 
diffusive scaling x ~ i( 1 + Q )/ Q . I n the limit a = 2, using 
2c 2 = a 2 we retrieve the previous result (TT5j) . while for 
a^Owe find regular Levy flight scaling x ~ t 1 /™ [f2l |. 




Figure 5: Time-averaged MSD as function of A for ten differ- 
ent trajectories of length T = 2000 with jump length correla- 
tion. We chose a 2 = 2. 



A. Weak ergodicity breaking 

The process with coupled jump lengths has the same 
waiting time for each jump. Could it still be subject to 
weak ergodicity breaking? Combining Eqs. (fTTj) and (fT6|) 
we obtain the equality 



A-U.Tj) = ^a 2 t- 



In the limit A « T wc find the scaling (^S 2 (A,T) 

A 2 T, contrasting the leading cubic behaviour (r 2 (t)) ~ t 3 
of the ensemble average. Thus, also the non-renewal 
CTRW with coupled jump lengths leads to a weak er- 
godicity breaking. Again, simulations corroborate the 
scaling with lag time A and overall process time T. In 
Fig. [5]we show the scatter between different trajectories. 



V. ANTICORRELATED WAITING TIMES 

Another way to introduce correlations is to consider 
diffusion on a network of nodes on each of which a differ- 
ent distribution is assumed. We illustrate this approach 
by anticorrelated waiting times. Namely we have two 
waiting time distributions, <J/i and ^2, such that $1 is 
centred around short waiting times and $2 around longer 
ones. We start by choosing a waiting time from one of 
the 'E'j and then pass from one waiting time distribution 
to the other according to the followings rules: (i) if we are 
in state ^1 and get a waiting time shorter than a preset 
t\, we shift to ^2 for the next step. Otherwise we remain 
with $1 . (ii) If we are in state ^2 and find a waiting time 
longer than a given ti, we shift to ^1, otherwise stay in 

*2- 

We can view this process as a diffusion in a network 
with two nodes. Being at node i at step n means that 
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the nth waiting time is extracted from distribution "I^. 
The transition probabilities between the two nodes are 
Pi = Sq ^i(t)dt and P 2 = /~ $ 2 (t)dt. The probability 
that the nth waiting time is chosen from distribution "fi 
is Pr(*i,n) = (l-Pi)Pr(*i,n-l)+P 2 Pr(* 2 ,n-l), 
and similarly for \E'2- These two equations can easily be 
solved. With a given initial condition we find 

Pr(^ ; n)= P2 + A "( PlPr ^' Q | J - P2Pr ^' » i (18) 

Pi + Pi 

with A = 1 — P\ — P 2 . For n — > oo we converge to the 
equilibrium probability Pr(*i) = P 2 /(Pl + Pi). If both 
distributions have a finite characteristic waiting time we 
see that at n -> oo, = (P 2 ^i) + Pi^ 2 ))/(Pi + P 2 )- 
Consider two limits: (i) if Pi = 1 and P 2 = 1 (i.e., 
<i — > oo and i 2 = 0) we change the waiting time at each 
step. Roughly, if we begin with the distribution 'J'i we ex- 
pect * = n ((V>i> + (V*» /2 + [1 - (-1)"] «^ 2 ) - (Vi» /4. 
If (i/ji) <C ("02), at the beginning we have a short wait- 
ing time alternating with a long one. After a long time 
t (tp2), however, we see a smooth curve for the mean 
squared displacement, see Fig. [51 (ii) If Pi <C 1, we have 
a classical CTRW governed by distribution tpi that after 
a while is somewhat modified by distribution tp 2 ■ At long 
times we find the same result as in case (i). This model 
can easily be extended to a wider network of nodes. 
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Figure 6: Anticorrelated CTRW for a 3D Gaussian walk. 
and ij)2 are normal distributions centred around and 100. 
The shift times are t\=ti = 50. 



Appendix A: Explicit derivation of the jump 
statistics 

We here show how the average number of jumps for 
correlated waiting time, and the position correlations are 
calculated. 



1. Waiting time correlation 



VI. CONCLUSIONS 

We demonstrated that the elusive coupled, non- 
renewal CTRW can indeed be applied to non-Brownian 
processes. Compared to previous models we believe that 
the idea of diffusion in the space of waiting times or jump 
length is quite intuitive and generic, such that this model 
will lend itself to a broad class of phenomena. In partic- 
ular, this approach allows us to consider Levy stable cor- 
relations in waiting times and jump lengths, significantly 
generalising previous results. We find that correlated 
waiting times lead to subdiffusion while correlations in 
the jump lengths give rise to superdiffusion, even when 
the waiting time or jump length increments are Gaus- 
sian. We also showed that anticorrelations in the long 
time limit produce normal diffusion. 

Both temporal and spatial correlations are demon- 
strated to lead to weak ergodicity breaking: the long 
time average of the mean squared displacement of a single 
trajectory shows different scaling behaviour than the cor- 
responding ensemble average. Surprisingly this is also 
the case for jump length correlations in which successive 
jumps are separated by constant waiting times. 
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To pass from Eq. (O to ([6|) we note that 
n— 1 n— 1 i n— 1 

E^ = EE^ = E^- ? )<^- ( A1 ) 

i=0 i=0 j=0 j=0 

For a stable law we know that the sum Y of independent, 
identically distributed random variables X fulfils [2l[ 

n / / n \ 1 / a \ 

Y = Y / a k X k ^Y(s)=X A - (A2) 



k=l 



with ELV( 



n—lf -\2 _ -2 



n — i 



see that in the Laplace domain 

n-l 



E"=i * = n ( n + !)(» + 2)/6 we 



n(n + l)(2n+ 1) 



k=0 



G 



s , (A3) 



leading to Eq. ([B]). The mean number of steps, (n(s)) 
Sn=o P( s ; n )i s then found to be 



1 °° I 
l (s)) =-X> 

71=1 \ 



i(n + l)(2n+ 1) 
6 



s . (A4) 



Due to its Gaussian nature we know that Sip(s) = 1 — ts+ 
o(s), where r is the mean waiting time for a Gaussian ran- 
dom variable of variance a. As we are only interested in 
the behaviour at long times, i.e., small s, we approximate 
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8ip(s) ~ exp(— ts) and derive the leading contribution of 
(n(s)) around s = 0. With n(n + l)(2n + l)/6 » n 3 /3, 

^ OO 

(„(«)) ~ _ £ exp (- [n 3 /3] 1/2 sr) . (A5) 



Approximating the sum by an integral, we finally obtain 
Eq. (0. For a Levy stable statistics of the increments ijji 
we use the characteristic function 5ip(uj) = exp (— c|w| a ) 
in Fourier space denoted by the frequency ui. We find 



i=0 V i=l / 

from which we obtain the average number of steps 
(nM) = f; n 1 = CXP H " + 1)Q|CU;|a) e- ^ 



n=0 



(A7) 

After approximation of the harmonic number X)i=i *° ~ 
n 1+Q /(l + a) and turning from sum to integral, we find 
after Fourier transform 



<n(t)> 



(o + l)V(«+i) /A a/(Q+1) 



2cos(a7r/[2(a + 1)]) \c / 
In the limit a — 2 we return to expression (J7J). 



(A8) 



Due to independence of the increments, (<5x(i) • 5x(j)) = 
^<r 2 5ij, we obtain the exact relation (fl~5l) . 

For the jump correlation (r(t) • r(t + r)) we start with 



(x(f) • x(i + r)) = ( x(t) • ( x(t) + ^ <5x(t + i) 



i=i 



(A10) 

As x(i) and <5x(t + i) are uncorrelated and of mean 0, we 
obtain (x(t) -x(t + r)) = |cr 2 <. This expression allows us 



(A6) 

to obtain the position correlation 



(r(t) -r(t+r)) = ( r(t) • r(t) + ]T x(i + i) ) . (All) 



Since x(i) have zero mean we obtain 



(r(f).r(f+r)) = ((r(i)) 2 )+^^(x(z). X (Aj)). (A12) 



2. Jump length correlation 



Assume the Gaussian distribution (|14|) for the jump 
displacement x(t) = r(t) — r(t — 1) in the fth jump, with 
initial condition x(0) = 0. The incremental change of the 
jump lengths is Sx(t) = x(t) — x(i — 1). Consequently, 
(x 2 (t)) = \a 2 t. We can then calculate the MSD (Ar(t) = 
r(t)-r(0)), 

= /(ec* +1 ) 

= E E(* + 1 - 0(* + 1 - i)<*x(0 ■ m?)> (A9) 
i=i j=i 



While the 



As (x(z) • x(f + j)) = |cr 2 i, we arrive at Eq. ([To]) . 

For a Levy distribution of the jump increments the 
position distribution of the overall process is stable with 
index a and scale factor c(t) - ci (1+Q)/a /(l + a) 1/a . If 
we concentrate on the x coordinate, the characteristic 

function becomes P(q, t) = exp qc(t) 

variance of this law diverges, one can calculate fractional 
moments of order < 5 < a 12], scaling like (|sg| (t)) ~ 
c(t) s . We therefore find the superdiffusive scaling x 2 ~ 
j2{i+ot)/a_ j n ^ e ij m j^ a — >. 2, this reproduces the cubic 
scaling x 2 ~ i 3 found for Gaussian correlations in the 
jump increments. For decreasing a the superdiffusion is 
enhanced, and in the limit of small a we find the scaling 
x 2 ~ t 2 / a of a regular Levy flight, i.e., a renewal process 
with a Levy stable jump length distribution A(£) of index 
a. 
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